
rm(list=ls())
gc()
setwd("/hpc/group/pbenfeylab/CheWei/")


as.numeric(system("awk '/MemFree/ {print $2}' /proc/meminfo", intern=TRUE))


suppressMessages(library(Seurat))
suppressMessages(library(tricycle))
suppressMessages(library(cowplot))
suppressMessages(library(ggplot2))
suppressMessages(library(scattermore))
suppressMessages(library(scater))
suppressMessages(library(cowplot))
suppressMessages(library(RColorBrewer))
suppressMessages(library(grid))
suppressMessages(library(gplots))
suppressMessages(library(tidyverse))

prep <- function(seu){
seu$time.anno.Li.crude <- gsub("_.*$","",seu$time.celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Distal Columella","Distal Columella_Columella",seu$time.celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Distal Lateral Root Cap","Distal Lateral Root Cap_Lateral Root Cap",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Proximal Columella","Proximal Columella_Columella",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Proximal Lateral Root Cap","Proximal Lateral Root Cap_Lateral Root Cap",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("^.*_","",seu$celltype.anno.Li.crude)
    return(seu)
}
plot_anno <- function(rc.integrated){
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno.Li.crude <- factor(rc.integrated$celltype.anno.Li.crude, levels = order[sort(match(unique(rc.integrated$celltype.anno.Li.crude),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno.Li.crude),order))]
p1 <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno.Li.crude", cols=color)
p2 <- DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno.Li.crude", order = c("Distal Columella","Proximal Columella","Distal Lateral Root Cap","Proximal Lateral Root Cap","Maturation","Elongation", "Transition Domain", "Proliferation Domain"),
        cols = c('#F7E7B0','#FFC400','#2B871F','#005E3B',"#deebf7", "#3182bd", '#fee0d2','#de2d26'))
p3 <- DimPlot(rc.integrated, reduction = "umap", group.by = "consensus.time.group", cols=brewer.pal(10,"Spectral"))
p4 <- DimPlot(rc.integrated, reduction = "umap", group.by = "branch.anno")
options(repr.plot.width=16, repr.plot.height=12)
gl <- lapply(list(p1, p2, p3, p4), ggplotGrob)
gwidth <- do.call(unit.pmax, lapply(gl, "[[", "widths"))
gl <- lapply(gl, "[[<-", "widths", value = gwidth)
gridExtra::grid.arrange(grobs=gl, ncol=2)
}


sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: AlmaLinux 9.3 (Shamrock Pampas Cat)

Matrix products: default
BLAS/LAPACK: /hpc/group/pbenfeylab/ch416/miniconda3/envs/seu4/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] forcats_0.5.2               stringr_1.5.1              
 [3] dplyr_1.1.3                 purrr_1.0.2                
 [5] readr_2.1.3                 tidyr_1.3.0                
 [7] tibble_3.2.1                tidyverse_1.3.2            
 [9] gplots_3.1.3                RColorBrewer_1.1-3         
[11] scater_1.26.1               scuttle_1.8.0              
[13] scattermore_1.2             ggplot2_3.4.4              
[15] cowplot_1.1.1               tricycle_1.6.0             
[17] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0
[19] Biobase_2.58.0              GenomicRanges_1.50.0       
[21] GenomeInfoDb_1.34.8         IRanges_2.32.0             
[23] S4Vectors_0.36.0            BiocGenerics_0.44.0        
[25] MatrixGenerics_1.10.0       matrixStats_1.1.0          
[27] SeuratObject_4.1.3          Seurat_4.1.1.9001          

loaded via a namespace (and not attached):
  [1] utf8_1.2.4                spatstat.explore_3.2-5   
  [3] reticulate_1.34.0         tidyselect_1.2.0         
  [5] RSQLite_2.3.1             AnnotationDbi_1.60.0     
  [7] htmlwidgets_1.6.2         BiocParallel_1.32.5      
  [9] Rtsne_0.16                munsell_0.5.0            
 [11] ScaledMatrix_1.6.0        codetools_0.2-19         
 [13] ica_1.0-3                 pbdZMQ_0.3-8             
 [15] future_1.33.0             miniUI_0.1.1.1           
 [17] withr_2.5.2               spatstat.random_3.2-1    
 [19] colorspace_2.1-0          progressr_0.14.0         
 [21] uuid_1.1-0                ROCR_1.0-11              
 [23] tensor_1.5                listenv_0.9.0            
 [25] repr_1.1.4                GenomeInfoDbData_1.2.9   
 [27] polyclip_1.10-6           bit64_4.0.5              
 [29] parallelly_1.36.0         vctrs_0.6.4              
 [31] generics_0.1.3            circular_0.4-95          
 [33] timechange_0.1.1          R6_2.5.1                 
 [35] ggbeeswarm_0.7.1          rsvd_1.0.5               
 [37] bitops_1.0-7              spatstat.utils_3.0-4     
 [39] cachem_1.0.8              DelayedArray_0.24.0      
 [41] assertthat_0.2.1          promises_1.2.1           
 [43] scales_1.2.1              googlesheets4_1.0.1      
 [45] beeswarm_0.4.0            gtable_0.3.4             
 [47] beachmat_2.14.0           globals_0.16.2           
 [49] goftest_1.2-3             rlang_1.1.2              
 [51] splines_4.2.2             lazyeval_0.2.2           
 [53] gargle_1.2.1              broom_1.0.2              
 [55] spatstat.geom_3.2-7       modelr_0.1.10            
 [57] reshape2_1.4.4            abind_1.4-5              
 [59] backports_1.4.1           httpuv_1.6.12            
 [61] tools_4.2.2               ellipsis_0.3.2           
 [63] ggridges_0.5.4            Rcpp_1.0.11              
 [65] plyr_1.8.9                base64enc_0.1-3          
 [67] sparseMatrixStats_1.10.0  zlibbioc_1.44.0          
 [69] RCurl_1.98-1.6            deldir_1.0-9             
 [71] pbapply_1.7-2             viridis_0.6.4            
 [73] zoo_1.8-12                haven_2.5.1              
 [75] ggrepel_0.9.4             cluster_2.1.4            
 [77] fs_1.6.3                  magrittr_2.0.3           
 [79] data.table_1.14.8         RSpectra_0.16-1          
 [81] reprex_2.0.2              lmtest_0.9-40            
 [83] RANN_2.6.1                googledrive_2.0.0        
 [85] mvtnorm_1.1-3             ggnewscale_0.4.8         
 [87] fitdistrplus_1.1-11       hms_1.1.2                
 [89] patchwork_1.1.3           mime_0.12                
 [91] evaluate_0.23             xtable_1.8-4             
 [93] readxl_1.4.1              fastDummies_1.7.3        
 [95] gridExtra_2.3             compiler_4.2.2           
 [97] KernSmooth_2.23-20        crayon_1.5.2             
 [99] htmltools_0.5.7           tzdb_0.3.0               
[101] later_1.3.1               lubridate_1.9.0          
[103] DBI_1.1.3                 dbplyr_2.2.1             
[105] MASS_7.3-58.3             boot_1.3-28.1            
[107] Matrix_1.5-4              cli_3.6.1                
[109] parallel_4.2.2            igraph_1.5.1             
[111] pkgconfig_2.0.3           sp_2.1-1                 
[113] IRdisplay_1.1             plotly_4.10.3            
[115] spatstat.sparse_3.0-3     xml2_1.3.3               
[117] vipor_0.4.5               XVector_0.38.0           
[119] rvest_1.0.3               digest_0.6.33            
[121] sctransform_0.4.1         RcppAnnoy_0.0.21         
[123] spatstat.data_3.0-3       Biostrings_2.66.0        
[125] cellranger_1.1.0          leiden_0.4.3             
[127] uwot_0.1.16               DelayedMatrixStats_1.20.0
[129] shiny_1.7.5.1             gtools_3.9.4             
[131] lifecycle_1.0.4           nlme_3.1-162             
[133] jsonlite_1.8.7            BiocNeighbors_1.16.0     
[135] viridisLite_0.4.2         fansi_1.0.5              
[137] pillar_1.9.0              lattice_0.21-8           
[139] KEGGREST_1.38.0           fastmap_1.1.1            
[141] httr_1.4.7                survival_3.4-0           
[143] glue_1.6.2                png_0.1-8                
[145] bit_4.0.5                 stringi_1.8.1            
[147] blob_1.2.4                RcppHNSW_0.5.0           
[149] BiocSingular_1.14.0       caTools_1.18.2           
[151] memoise_2.0.1             IRkernel_1.3.1.9000      
[153] irlba_2.3.5.1             future.apply_1.11.0


cc.genes.goldy <- read.csv("./tricycle/Goldy_2021_CycB1-1-sorted_Enriched_Depleted_G2M_rm-proto-genes.csv")
cc.genes.zhang <- read.csv("./tricycle/Zhang_2021_Dev_cell_Core_cell_cycle_genes.csv")
cc.genes.amigo.go <- as.character(read.table("./tricycle/GO_cell_cycle_annotation_arabidopsis.txt", header=F)$V1)
cc.genes.f <- unique(c(cc.genes.zhang$Gene, cc.genes.goldy$gene, cc.genes.amigo.go))
length(cc.genes.f)
pp.genes <- as.character(read.table("./CW_data/escoring/Root_sc/Protoplasting_DEgene_FC2_list.txt", header=F)$V1)


table(cc.genes.zhang$Phase)

G1/G0    G2     M     S 
   18    95    77    53


head(cc.genes.zhang)


## Cell Cyle genes used
cc.gene.list <- as.list(cc.genes.f)


rc.integrated <- readRDS("./scRNA-seq/Integrated_Objects/rc.integrated_18S_WT_cell_cycle_seu4_245genes_Goldy_Zhang_AMIGO_20230825.rds")


rownames(rc.integrated)


length(rownames(rc.integrated))


write.csv(data.frame("GeneID"=rownames(rc.integrated)),"./tradeseq/245_cell_cycle_related_genes_for_reference.csv")


## Prepare New Ref (All cells)
corrected.m <- rc.integrated@assays$integrated@data


# run PCA on the batch effects corrected matrix and get the rotaions scores for the top 2 PCs
pca.m <- scater::calculatePCA(corrected.m, ntop = 500)


new.ref <- attr(pca.m, 'rotation')[, seq_len(2)]


saveRDS(new.ref, "./tradeseq/Tricycle_Reference_Arabidopsis_Root.rds")


unique(rc.integrated$orig.ident)


#suppressMessages(rc.integrated <- RunUMAP(rc.integrated, reduction = "pca", dims = 1:8))


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "seurat_clusters", label = TRUE)+theme(legend.position="none")

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "seurat_clusters", label = TRUE)+theme(legend.position="none")

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "orig.ident", label = FALSE, pt.size=1, order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "orig.ident", label = FALSE, pt.size=1, order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Xylem Pole Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno.Li <- factor(rc.integrated$celltype.anno.Li, levels = order[sort(match(unique(rc.integrated$celltype.anno.Li),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno.Li),order))]


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno.Li", label = FALSE, cols=color, pt.size=1, order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "celltype.anno.Li", label = FALSE, cols=color, pt.size=1, order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "time.anno.Li", label = FALSE, pt.size=1, order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


options(repr.plot.width=8, repr.plot.height=8)
DimPlot(rc.integrated, reduction = "pca", group.by = "consensus.time.group", label = FALSE, pt.size=1, order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


## S-phase cluster
sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==0)]$celltype.anno.Li))

     Quiescent Center                Phloem Phloem Pole Pericycle 
                   38                   275                   323 
            Columella                 Xylem                Cortex 
                  331                   795                  1106 
           Procambium  Xylem Pole Pericycle            Endodermis 
                 1426                  1439                  1592 
         Atrichoblast           Trichoblast      Lateral Root Cap 
                 1646                  2546                  2879


## S-phase cluster
sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==0)]$time.anno.Li))

  Distal Lateral Root Cap          Distal Columella        Proximal Columella 
                      136                       161                       170 
               Maturation         Transition Domain Proximal Lateral Root Cap 
                      763                      1900                      2743 
     Proliferation Domain                Elongation 
                     3610                      4913


## G2M-phase cluster
rev(sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==4)]$celltype.anno.Li)))

     Lateral Root Cap          Atrichoblast            Endodermis 
                 2981                  1013                   790 
          Trichoblast                Cortex                 Xylem 
                  550                   429                   284 
            Columella  Xylem Pole Pericycle            Procambium 
                  275                   252                   249 
               Phloem Phloem Pole Pericycle      Quiescent Center 
                  233                   100                     8


## G2M-phase cluster
rev(sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==4)]$time.anno.Li)))

Proximal Lateral Root Cap      Proliferation Domain                Elongation 
                     2863                      1840                      1324 
               Maturation         Transition Domain          Distal Columella 
                      533                       211                       185 
  Distal Lateral Root Cap        Proximal Columella 
                      118                        90


## G2M-phase cluster
rev(sort(table(rc.integrated[,which(rc.integrated$seurat_clusters==4)]$time.celltype.anno.Li)))

                       Proximal Lateral Root Cap 
                                            2863 
               Proliferation Domain_Atrichoblast 
                                             828 
                 Proliferation Domain_Endodermis 
                                             509 
                     Proliferation Domain_Cortex 
                                             264 
                Proliferation Domain_Trichoblast 
                                             210 
                 Elongation_Xylem Pole Pericycle 
                                             196 
                           Elongation_Endodermis 
                                             195 
                                Distal Columella 
                                             185 
                           Elongation_Procambium 
                                             176 
          Elongation_Metaphloem & Companion Cell 
                                             153 
                   Transition Domain_Trichoblast 
                                             144 
                           Elongation_Protoxylem 
                                             128 
                           Maturation_Protoxylem 
                                             123 
                               Elongation_Cortex 
                                             118 
                         Distal Lateral Root Cap 
                                             118 
                          Elongation_Trichoblast 
                                             114 
                         Elongation_Atrichoblast 
                                              99 
                              Proximal Columella 
                                              90 
                Elongation_Phloem Pole Pericycle 
                                              90 
                           Maturation_Endodermis 
                                              85 
                          Maturation_Trichoblast 
                                              82 
                           Maturation_Procambium 
                                              71 
                         Maturation_Atrichoblast 
                                              57 
                 Maturation_Xylem Pole Pericycle 
                                              41 
                               Maturation_Cortex 
                                              33 
                  Transition Domain_Atrichoblast 
                                              29 
                          Elongation_Protophloem 
                                              29 
          Maturation_Metaphloem & Companion Cell 
                                              28 
                            Elongation_Metaxylem 
                                              26 
   Transition Domain_Metaphloem & Companion Cell 
                                              14 
                        Transition Domain_Cortex 
                                              14 
       Proliferation Domain_Xylem Pole Pericycle 
                                              12 
           Proliferation Domain_Quiescent Center 
                                               8 
                Maturation_Phloem Pole Pericycle 
                                               8 
                  Proliferation Domain_Metaxylem 
                                               5 
                          Maturation_Protophloem 
                                               5 
          Transition Domain_Xylem Pole Pericycle 
                                               3 
                   Transition Domain_Protophloem 
                                               2 
         Transition Domain_Phloem Pole Pericycle 
                                               2 
Proliferation Domain_Metaphloem & Companion Cell 
                                               2 
                    Transition Domain_Procambium 
                                               1 
                     Transition Domain_Metaxylem 
                                               1 
                    Transition Domain_Endodermis 
                                               1 
                 Proliferation Domain_Protoxylem 
                                               1 
                 Proliferation Domain_Procambium 
                                               1


pc1r <- as.data.frame(rc.integrated[['pca']][,1]) %>% arrange(desc(PC_1)) %>% rownames(.)


pc1r[1:9]


## PC1 positive
FeaturePlot(rc.integrated, reduction="pca", features = c(pc1r[1:9]), cols = c("#EEEEEE","#FF4040"), order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


## PC1 positive
FeaturePlot(rc.integrated, reduction="umap", features = c(pc1r[1:9]), cols = c("#EEEEEE","#FF4040"), order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


pc1r[1:20]


## PC1 negative
FeaturePlot(rc.integrated, reduction="pca", features = c(rev(pc1r)[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


pc2r <- as.data.frame(rc.integrated[['pca']][,2]) %>% arrange(desc(PC_2)) %>% rownames(.)


pc2r[1:20]


rev(pc2r)[1:20]


## PC2 positive
FeaturePlot(rc.integrated, reduction="pca", features = c(pc2r[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


## PC2 negative
FeaturePlot(rc.integrated, reduction="pca", features = c(rev(pc2r)[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


## PC2 negative
FeaturePlot(rc.integrated, reduction="umap", features = c(rev(pc2r)[1:9]), cols = c("#EEEEEE","#FF4040"),order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


pc3r <- as.data.frame(rc.integrated[['pca']][,3]) %>% arrange(desc(PC_3)) %>% rownames(.)


pc3r[1:20]


rev(pc3r)[1:20]


zscore <- function(x){(x-mean(x))/sd(x)}


cc.genes.zhang <- cc.genes.zhang[match(intersect(rownames(rc.integrated),cc.genes.zhang$Gene),cc.genes.zhang$Gene),] %>% arrange(Phase)


head(cc.genes.zhang)


head(cc.genes.zhang %>% filter(Phase == "S"))


msc <- c()
suppressWarnings(
for (i in as.character(unique(cc.genes.zhang$Phase))){
    if (length(cc.genes.zhang[which(cc.genes.zhang$Phase== i),]$Gene)>1){
    msc <- cbind(msc, as.numeric(apply(apply(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.zhang[which(cc.genes.zhang$Phase== i),]$Gene,], 1, zscore), 1, mean)))       
    } else {
    msc <- cbind(msc, as.numeric(zscore(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.zhang[which(cc.genes.zhang$Phase== i),]$Gene,])))      
    }

}
)
colnames(msc) <- as.character(unique(cc.genes.zhang$Phase))
rownames(msc) <- colnames(rc.integrated)


head(msc)


#rc.integrated$G1_sig <- msc[,1]
rc.integrated$G2_sig <- msc[,1]
rc.integrated$M_sig <- msc[,2]
rc.integrated$S_sig <- msc[,3]


#FeaturePlot(rc.integrated, reduction="umap", features = "G1_sig", pt.size=1,order=TRUE)


FeaturePlot(rc.integrated, reduction="umap", features = "S_sig", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


FeaturePlot(rc.integrated, reduction="umap", features = "G2_sig", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


FeaturePlot(rc.integrated, reduction="umap", features = "M_sig", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


#FeaturePlot(rc.integrated, reduction="pca", features = "G1_sig", pt.size=1,order=TRUE)


FeaturePlot(rc.integrated, reduction="pca", features = "S_sig", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


FeaturePlot(rc.integrated, reduction="pca", features = "G2_sig", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


FeaturePlot(rc.integrated, reduction="pca", features = "M_sig", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


head(cc.genes.goldy)


cc.genes.goldy <- cc.genes.goldy[match(intersect(rownames(rc.integrated),cc.genes.goldy$gene),cc.genes.goldy$gene),] %>% arrange(phase)


msc <- c()
suppressWarnings(
for (i in as.character(unique(cc.genes.goldy$phase))){
    if (length(cc.genes.goldy[which(cc.genes.goldy$phase== i),]$gene)>1){
    msc <- cbind(msc, as.numeric(apply(apply(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.goldy[which(cc.genes.goldy$phase== i),]$gene,], 1, zscore), 1, mean)))       
    } else {
    msc <- cbind(msc, as.numeric(zscore(as.matrix(rc.integrated@assays$integrated@data)[cc.genes.goldy[which(cc.genes.goldy$phase== i),]$gene,])))      
    }

}
)
colnames(msc) <- as.character(unique(cc.genes.goldy$phase))
rownames(msc) <- colnames(rc.integrated)


head(msc)


rc.integrated$Depleted_G2M <- msc[,1]
rc.integrated$Enriched_G2M <- msc[,2]


FeaturePlot(rc.integrated, reduction="umap", features = "Depleted_G2M", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


FeaturePlot(rc.integrated, reduction="umap", features = "Enriched_G2M", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


FeaturePlot(rc.integrated, reduction="pca", features = "Depleted_G2M", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


FeaturePlot(rc.integrated, reduction="pca", features = "Enriched_G2M", pt.size=1,order=TRUE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


rc.integrated$G2M_clust <- rep("others", ncol(rc.integrated))
rc.integrated$G2M_clust[which(rc.integrated$seurat_clusters==4)] = "G2M"


DimPlot(rc.integrated, reduction="umap", group = "G2M_clust", pt.size=1,order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


DimPlot(rc.integrated, reduction="pca", group = "G2M_clust", pt.size=1,order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


rc.integrated$S_clust <- rep("others", ncol(rc.integrated))
rc.integrated$S_clust[which(rc.integrated$seurat_clusters==0)] = "S"


DimPlot(rc.integrated, reduction="umap", group = "S_clust", pt.size=1,order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


DimPlot(rc.integrated, reduction="pca", group = "S_clust", pt.size=1,order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


new.ref <- readRDS("./tradeseq/Tricycle_Reference_Arabidopsis_Root.rds")


sce <- as.SingleCellExperiment(rc.integrated)


sce$geno <- rep("WT",length(sce$orig.ident))


## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)

sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi

## Plot
options(repr.plot.width=9, repr.plot.height=4)
print(plot_ccposition_den(sce$tricyclePosition,
                sce$geno, 'genotype',
                bw = 10, fig.title = "Kernel density of \u03b8", palette.v=colorRampPalette(brewer.pal(11, "Spectral"))(18)) +
  theme_bw(base_size = 14)+ geom_vline(xintercept = 1*pi, linetype="dotted", 
                color = "red", size=1.5)+ geom_vline(xintercept = 1.55*pi, linetype="dotted", 
                color = "red", size=1.5))

The number of projection genes found in the new data is 245.

Warning message:
“Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.”


grDevices::cairo_pdf("./pdfs/Density_Plot_CC.pdf",width=9, height=4)
print(plot_ccposition_den(sce$tricyclePosition,
                sce$geno, 'genotype',
                bw = 10, fig.title = "Kernel density of \u03b8", palette.v=colorRampPalette(brewer.pal(11, "Spectral"))(18)) +
  theme_bw(base_size = 14)+ geom_vline(xintercept = 1*pi, linetype="dotted", 
                color = "red", size=1.5)+ geom_vline(xintercept = 1.55*pi, linetype="dotted", 
                color = "red", size=1.5))
dev.off()


p <- plot_emb_circle_scale(sce, dimred = 1,
                           point.size = 5, point.alpha = 0.9) +
  theme_bw(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
options(repr.plot.width=8, repr.plot.height=4)
plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.6))


grDevices::cairo_pdf("./pdfs/PCA_CC.pdf",width=8, height=4)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.6)))
dev.off()


## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)

sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
sce$tricyclePosition <- sce$tricyclePosition*pi

## Plot
options(repr.plot.width=9, repr.plot.height=9)
print(plot_ccposition_den(sce$tricyclePosition,
                sce$orig.ident, 'sample',
                bw = 10, fig.title = "Kernel density of \u03b8", palette.v=colorRampPalette(brewer.pal(11, "Spectral"))(18)) +
  theme_bw(base_size = 14)+ geom_vline(xintercept = 1*pi, linetype="dotted", 
                color = "red", size=1.5)+ geom_vline(xintercept = 1.55*pi, linetype="dotted", 
                color = "red", size=1.5))

The number of projection genes found in the new data is 245.


## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    rc.integrated$tricyclePosition <- sce$tricyclePosition
    rc.integrated$tricycleCCStage <- tricycleCCStage


saveRDS(rc.integrated,"./scRNA-seq/Integrated_Objects/rc.integrated_18S_WT_cell_cycle_seu4_245genes_Goldy_Zhang_AMIGO_20230825.rds")


## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 1, point.alpha = 0.9) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))


DimPlot(rc.integrated, reduction="umap", group = "tricycleCCStage", pt.size=1,order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


DimPlot(rc.integrated, reduction="pca", group = "tricycleCCStage", pt.size=1,order=FALSE)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`


rc.integrated$CCclust <- rep("others", ncol(rc.integrated))
rc.integrated$CCclust[which(rc.integrated$seurat_clusters==0)] = "S"
rc.integrated$CCclust[which(rc.integrated$seurat_clusters==4)] = "G2M"
table(rc.integrated$tricycleCCStage, rc.integrated$CCclust)

       
           G2M others      S
  G1/G0    920 154117   1510
  G2M     6101   4126    117
  S        143  20007  12769


scater::plotReducedDim(sce, dimred = "tricycleEmbedding") +
    labs(x = "Projected PC1", y = "Projected PC2") +
    ggtitle(sprintf("Projected cell cycle space (n=%d)",
                    ncol(sce))) +
    theme_bw(base_size = 14)


p <- plot_emb_circle_scale(sce, dimred = "tricycleEmbedding",
                           point.size = 3.5, point.alpha = 0.9) +
  theme_bw(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.4))


## Total RNA content
#fit.l <- fit_periodic_loess(sce$tricyclePosition,
#                            log2(sce$nCount_RNA),
#                            plot = TRUE, span=0.3,
#                       x_lab = "Cell cycle position \u03b8", y_lab = "log2(Total UMIs)",
#                       fig.title = paste0("Expression of key gene along \u03b8 (n=",
#                                          ncol(sce), ")"))

#fit <- fit.l$pred.df
#options(repr.plot.width=8, repr.plot.height=6)
#ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "log2(total UMIs)", title = "") + 
#            scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0, 
#                pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0, 
#                0.5, 1, 1.5, 2), "π"), name = "θ")+
#  theme_bw(base_size = 14)+
#  theme(legend.title=element_blank())


table(rc.integrated$time.anno.Li)

         Distal Columella   Distal Lateral Root Cap                Elongation 
                     6526                      4904                     83690 
               Maturation      Proliferation Domain        Proximal Columella 
                    21054                     21913                      3063 
Proximal Lateral Root Cap         Transition Domain 
                    27495                     31165


sub <- subset(rc.integrated, cells=colnames(rc.integrated)[which(rc.integrated$time.anno.Li=="Proliferation Domain"|
                                                                rc.integrated$time.anno.Li=="Proximal Columella"|
                                                                rc.integrated$time.anno.Li=="Proximal Lateral Root Cap")])

sub

An object of class Seurat 
57118 features across 52471 samples within 3 assays 
Active assay: integrated (245 features, 245 variable features)
 2 layers present: data, scale.data
 2 other assays present: RNA, SCT
 3 dimensional reductions calculated: pca, umap, umap_2D


table(sub$time.anno.Li)

     Proliferation Domain        Proximal Columella Proximal Lateral Root Cap 
                    21913                      3063                     27495


# Run PCA
sub <- RunPCA(sub, npcs = 10, verbose = FALSE, approx = TRUE)

# Find nearest neighbors
#rc.integrated <- FindNeighbors(rc.integrated, reduction = "pca",dims = 1:50)

Warning message:
“Number of dimensions changing from 50 to 10”
Warning message:
“The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session”
20:27:45 UMAP embedding parameters a = 0.9922 b = 1.112

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

20:27:45 Read 52471 rows and found 5 numeric columns

20:27:45 Using Annoy for neighbor search, n_neighbors = 30

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

20:27:45 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

20:27:57 Writing NN index file to temp file /tmp/Rtmpu0VjyG/file317db91a27d9c1

20:27:57 Searching Annoy index using 1 thread, search_k = 3000

20:28:40 Annoy recall = 100%

20:28:42 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

20:28:44 5 smooth knn distance failures

20:28:51 Initializing from normalized Laplacian + noise (using RSpectra)

20:28:59 Commencing optimization for 200 epochs, with 1889004 positive edges

20:29:42 Optimization finished

Warning message:
“Key ‘UMAP_’ taken, using ‘umap_’ instead”


# Run 2D UMAP
sub <- RunUMAP(sub, reduction = "pca", dims = 1:3, n.neighbors = 30, min.dist=0.3)

20:40:15 UMAP embedding parameters a = 0.9922 b = 1.112

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

20:40:15 Read 52471 rows and found 3 numeric columns

20:40:15 Using Annoy for neighbor search, n_neighbors = 30

Found more than one class "dist" in cache; using the first, from namespace 'spam'

Also defined by ‘BiocGenerics’

20:40:15 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

20:40:28 Writing NN index file to temp file /tmp/Rtmpu0VjyG/file317db92623cbe

20:40:28 Searching Annoy index using 1 thread, search_k = 3000

20:40:57 Annoy recall = 100%

20:40:59 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

20:41:00 18 smooth knn distance failures

20:41:06 Initializing from normalized Laplacian + noise (using RSpectra)

20:41:40 Commencing optimization for 200 epochs, with 1402752 positive edges

20:42:13 Optimization finished

Warning message:
“Key ‘UMAP_’ taken, using ‘umap_’ instead”


my_cols = brewer.pal(8,"Set2")


options(repr.plot.width=10, repr.plot.height=8)
DimPlot(sub, reduction = "umap", group.by = "time.anno.Li",cols=alpha(my_cols,0.8), label = FALSE, pt.size=1)


sce <- as.SingleCellExperiment(sub)


# dims = 3
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 2, point.alpha = 0.8) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))


# dims = 4
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 2, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))


# dims = 5
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 2, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))


# dims = 5
options(repr.plot.width=8, repr.plot.height=8)
## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 2, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))


run_tricycle <- function(sample){
    seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
    sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)
    
    sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
    
    ## Write annotations as csv
    write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
    #saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}


sample.list <- c("sc_130","sc_131","sc_132","sc_133","sc_134","sc_135","sc_136","sc_137")


for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}

The number of projection genes found in the new data is 230.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 237.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


sample.list <- c("sc_5")


for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


sample.list <- c("sc_1","sc_2","sc_111","sc_112","sc_113","sc_114")


for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 243.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


sample.list <- c("sc_151","sc_152","sc_153","sc_154")


for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


sample.list <- c("pp1","col0","tnw2","sc_1","sc_9_at","sc_30", "sc_31", "sc_37", "sc_40", "sc_51")


for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 243.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 244.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 241.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


run_tricycle <- function(sample){
    seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
    sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)
    
    sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.2
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- (sce$tricyclePosition-2)*-1
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
    
    ## Write annotations as csv
    write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
    #saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}


sample.list <- c("dc1","dc2","tnw1","briT","briTR","sc_10_at","sc_11","sc_12","sc_20","sc_43","sc_44","sc_45","sc_46","sc_47","sc_48","sc_49","sc_50","sc_71",
                 "sc_122","sc_123","sc_124","sc_125","sc_126","sc_127","sc_128","sc_129","sc_130","sc_131","sc_132","sc_133","sc_134","sc_135","sc_136","sc_137",
                 "sc_155","sc_156","sc_157","sc_158","sc_159","sc_160","sc_161","sc_162","sc_163","sc_164","sc_165","sc_166","sc_167","sc_168","sc_169",
                 "sc_170","sc_171","sc_172","sc_173","sc_175","sc_176","sc_177","sc_178","sc_179","sc_180","sc_186","sc_187","sc_189","sc_190","sc_191",
                 "sc_215","sc_216","sc_217","sc_218","sc_220","sc_221","sc_222","sc_223","sc_224","sc_225","sc_226","sc_227","sc_228","sc_229","sc_235")


for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/","dc1","_COPILOT.rds"))


sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)

The number of projection genes found in the new data is 245.


sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


run_tricycle <- function(sample){
    seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
    sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)
    
    sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
    
    ## Write annotations as csv
    write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
    #saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}


mtx <- read.table('./scRNA-seq/Published_Data/GSE140778_20190115_data.txt')


meta <- read.csv('./scRNA-seq/Published_Data/GSE140778_20190115_meta_data.csv')


sum(!is.na(meta$Cluster_PSE))


head(mtx)


rownames(meta) <- meta$ID


seu <- CreateSeuratObject(mtx, meta.data=meta)

seu

An object of class Seurat 
32824 features across 1375 samples within 1 assay 
Active assay: RNA (32824 features, 0 variable features)


head(seu@meta.data)


seu <- NormalizeData(seu)
seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = nrow(seu))
all.genes <- rownames(seu)
seu <- ScaleData(seu, features = all.genes)
seu <- RunPCA(seu,features = VariableFeatures(object = seu))

Centering and scaling data matrix

PC_ 1 
Positive:  AT5G03240, AT3G53990, AT4G05320, AT1G70670, AT4G05050, AT1G75280, AT4G39090, AT3G22400, AT4G21980, AT1G29520 
	   AT3G03200, AT1G77700, AT5G48060, AT1G54400, AT3G06420, AT3G58720, AT4G33490, AT2G46630, AT1G29380, AT1G08160 
	   AT1G79480, AT5G07220, AT5G49320, AT5G18970, AT5G03510, AT2G41200, AT2G42820, AT1G76970, AT1G63310, AT1G61760 
Negative:  AT3G04400, AT3G24830, AT3G25520, AT5G20290, AT5G02960, AT1G08360, AT5G16130, AT1G18540, AT3G02080, AT3G49910 
	   AT4G27090, AT2G34480, AT5G15200, AT2G47610, AT5G10360, AT4G00100, AT4G34670, AT5G23740, AT2G39460, AT3G04920 
	   AT3G53020, AT1G69620, AT2G18020, AT2G27530, AT1G67430, AT3G14600, AT3G11940, AT3G09200, AT2G19730, AT4G18100 
PC_ 2 
Positive:  AT5G42050, AT1G20450, AT1G19180, AT1G20440, AT2G47730, AT1G78380, AT4G30960, AT3G10985, AT5G10830, AT5G11670 
	   AT5G45110, AT1G78660, AT5G63790, AT1G04250, AT2G29440, AT5G08790, AT4G16380, AT4G34150, AT5G64260, AT5G19110 
	   AT2G24850, AT4G35090, AT4G19200, AT3G52470, AT5G06320, AT4G21990, AT1G55920, AT1G78080, AT1G76680, AT4G36040 
Negative:  AT5G62960, AT1G80940, AT2G38550, AT1G05610, AT2G22730, AT1G31880, AT5G03530, AT1G55760, AT2G18380, AT1G62045 
	   AT1G79430, AT5G40580, AT1G63310, AT1G68000, AT1G61760, AT3G49360, AT4G23900, AT1G69970, AT2G20515, AT2G36400 
	   AT4G17510, AT5G26940, AT1G14900, AT5G04060, AT1G06490, AT2G32280, AT5G17260, AT5G48060, AT4G24250, AT1G79480 
PC_ 3 
Positive:  AT1G03220, AT1G48920, AT1G66580, AT5G52470, AT3G15000, AT5G20160, AT3G10690, AT1G74560, AT4G12600, AT1G56110 
	   AT5G27120, AT4G26110, AT2G38810, AT5G25370, AT4G25630, AT1G68230, AT3G05060, AT3G07770, AT3G44750, AT5G22650 
	   AT3G45970, AT1G52930, AT3G23990, AT3G58350, AT3G56070, AT2G19540, AT3G57150, AT4G15910, AT3G03950, AT1G26680 
Negative:  AT1G78040, AT5G60920, AT5G59880, AT2G45470, AT3G12610, AT5G39320, AT5G05170, AT5G49720, AT1G53500, AT1G75680 
	   AT1G08200, AT4G12730, AT5G12250, AT2G20750, AT1G26850, AT1G70490, AT1G19835, AT4G31590, AT3G12710, AT1G10630 
	   AT1G76160, AT4G34490, AT1G45688, AT4G32410, AT1G76670, AT5G59290, AT5G64740, AT5G47770, AT5G25100, AT3G07010 
PC_ 4 
Positive:  AT1G10760, AT4G01610, AT1G54030, AT1G32860, AT2G34020, AT1G18210, AT2G29980, AT3G55760, AT5G48300, AT3G12490 
	   AT4G09020, AT4G31290, AT1G21340, AT3G08690, AT1G56680, AT2G15490, AT2G32430, AT2G04025, AT4G00490, AT5G08590 
	   AT2G28900, AT3G46440, AT5G18470, AT2G41800, AT4G34970, AT1G72490, AT1G13620, AT4G23060, AT3G18230, AT1G19850 
Negative:  AT4G14130, AT1G62480, AT4G21960, AT3G53420, AT2G02130, AT5G59090, AT1G52760, AT3G21770, AT2G43910, AT3G04730 
	   AT4G29100, AT3G19450, AT3G23820, AT3G51950, AT5G56540, AT3G14310, AT4G11210, AT3G56240, AT2G30520, AT3G16330 
	   AT5G24800, AT4G19410, AT2G31083, AT5G51780, AT2G44790, AT1G10682, AT3G16340, AT3G11930, AT5G25890, AT4G36710 
PC_ 5 
Positive:  AT4G15960, AT1G70670, AT3G02470, AT4G18710, AT2G32150, AT2G39420, AT2G36650, AT3G14780, AT4G37730, AT5G16650 
	   AT3G47160, AT4G35450, AT4G35790, AT2G38800, AT5G13740, AT3G11280, AT4G39090, AT2G19110, AT5G64572, AT1G75380 
	   AT5G49320, AT1G80040, AT1G44170, AT5G18980, AT3G60680, AT3G63000, AT3G14350, AT4G21110, AT2G39080, AT2G21270 
Negative:  ATMG01390, ATMG00020, AT5G18610, AT2G07698, AT1G77690, AT5G23860, ATMG01190, AT1G04680, AT5G26270, AT3G02230 
	   AT2G04780, ATMG01360, AT1G64640, AT3G54810, AT1G69690, AT2G14890, ATCG00500, AT2G20750, ATMG00660, AT2G47650 
	   ATCG01180, ATCG00470, AT3G48185, AT2G21160, AT3G58120, ATCG00950, AT5G15650, ATMG00090, AT2G34300, AT1G21750


seu <- FindNeighbors(seu, dims = 1:10)
seu <- FindClusters(seu, resolution = 0.5)

Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 1375
Number of edges: 33789

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9051
Number of communities: 12
Elapsed time: 0 seconds


seu <- RunUMAP(seu, dims = 1:10)

Warning message:
“The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session”
10:49:41 UMAP embedding parameters a = 0.9922 b = 1.112

10:49:41 Read 1375 rows and found 10 numeric columns

10:49:41 Using Annoy for neighbor search, n_neighbors = 30

10:49:41 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

10:49:41 Writing NN index file to temp file /tmp/RtmpSzDebc/file16ca8b2394d636

10:49:41 Searching Annoy index using 1 thread, search_k = 3000

10:49:42 Annoy recall = 100%

10:49:43 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

10:49:46 Initializing from normalized Laplacian + noise (using RSpectra)

10:49:46 Commencing optimization for 500 epochs, with 49356 positive edges

10:49:49 Optimization finished


# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(seu, reduction = "umap")


head(seu@meta.data)

seu

An object of class Seurat 
32824 features across 1375 samples within 1 assay 
Active assay: RNA (32824 features, 27218 variable features)
 2 dimensional reductions calculated: pca, umap


#seu[["umap"]]@cell.embeddings[,1] <- seu[["umap"]]@cell.embeddings[,1]*-1
seu[["umap"]]@cell.embeddings[,2] <- seu[["umap"]]@cell.embeddings[,2]*-1
#u2 <- rc.integrated@reductions$umap@cell.embeddings[,1]
#u1 <- rc.integrated@reductions$umap@cell.embeddings[,2]
#rc.integrated@reductions$umap@cell.embeddings[,1] <- u1
#rc.integrated@reductions$umap@cell.embeddings[,2] <- u2


DimPlot(seu, reduction = "umap", group.by = "Cluster_PSE")


FeaturePlot(seu, reduction = "umap", features = "Pseudotime_PSE")


sce <- as.SingleCellExperiment(seu)

## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)

The number of projection genes found in the new data is 245.


sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 4, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))


DimPlot(seu, reduction = "umap", group.by = "tricycleCCStage")


seu$Pseudotime_PSE


seu$Pseudotime_PSE_rank <- rank(seu$Pseudotime_PSE)


FeaturePlot(seu, reduction = "umap", features = "Pseudotime_PSE_rank")


saveRDS(seu,"./tricycle/seu_Phloem_1375_cells.rds")


seu_PSE <- subset(seu, cells=colnames(seu)[!is.na(seu$Cluster_PSE)])


seu_PSE

An object of class Seurat 
32824 features across 758 samples within 1 assay 
Active assay: RNA (32824 features, 27218 variable features)
 2 dimensional reductions calculated: pca, umap


DimPlot(seu_PSE, reduction = "umap", group.by = "tricycleCCStage")


FeaturePlot(seu_PSE, reduction = "umap", features = "Pseudotime_PSE")


FeaturePlot(seu_PSE, reduction = "umap", features = "Pseudotime_PSE_rank")


saveRDS(seu_PSE,"./tricycle/seu_PSE_758_cells.rds")


## All cells (PSE, MSE, Procambium)
pltsg <- function(gene){
    keygene.idx <- which(rownames(seu@assays$RNA@data) == gene)
    fit.l <- fit_periodic_loess(seu$tricyclePosition*pi,
                            seu@assays$RNA@data[keygene.idx,],
                            plot = FALSE, span=0.3,
                       x_lab = "Cell cycle position \u03b8", y_lab = "log2(keygene)")
    fit <- fit.l$pred.df
    options(repr.plot.width=8, repr.plot.height=6)
    print(ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "smoothed gene expression", title = gene) + 
                scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0, 
                    pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0, 
                    0.5, 1, 1.5, 2), "π"), name = "θ")+
      #              scale_x_continuous(limits = c(0, 200), breaks = c(0, 
      #              50, 100, 150, 200), labels = paste0(c(0, 
      #              0.5, 1, 1.5, 2), "π"), name = "θ")+
      theme_bw(base_size = 14)+
      theme(legend.title=element_blank()))
}


## DWF4
pltsg('AT3G50660')


## OPS
pltsg('AT3G09070')


## CPD
pltsg('AT5G05690')


## Only PSE
pltsg <- function(gene){
    keygene.idx <- which(rownames(seu_PSE@assays$RNA@data) == gene)
    fit.l <- fit_periodic_loess(seu_PSE$tricyclePosition*pi,
                            seu_PSE@assays$RNA@data[keygene.idx,],
                            plot = FALSE, span=0.3,
                       x_lab = "Cell cycle position \u03b8", y_lab = "log2(keygene)")
    fit <- fit.l$pred.df
    options(repr.plot.width=8, repr.plot.height=6)
    print(ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "smoothed gene expression", title = gene) + 
                scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0, 
                    pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0, 
                    0.5, 1, 1.5, 2), "π"), name = "θ")+
      #              scale_x_continuous(limits = c(0, 200), breaks = c(0, 
      #              50, 100, 150, 200), labels = paste0(c(0, 
      #              0.5, 1, 1.5, 2), "π"), name = "θ")+
      theme_bw(base_size = 14)+
      theme(legend.title=element_blank()))
}


## DWF4
pltsg('AT3G50660')


## OPS
pltsg('AT3G09070')


## CPD
pltsg('AT5G05690')


seu_mer_PSE <- subset(seu_PSE, cells=colnames(seu_PSE)[which(seu_PSE$Pseudotime_PSE_rank <= 410)])


seu_mer_PSE

An object of class Seurat 
32824 features across 410 samples within 1 assay 
Active assay: RNA (32824 features, 27218 variable features)
 2 dimensional reductions calculated: pca, umap


saveRDS(seu_mer_PSE,"./tricycle/seu_mer_PSE_410_cells.rds")


## Only meristem PSE
pltsg <- function(gene){
    keygene.idx <- which(rownames(seu_mer_PSE@assays$RNA@data) == gene)
    fit.l <- fit_periodic_loess(seu_mer_PSE$tricyclePosition*pi,
                            seu_mer_PSE@assays$RNA@data[keygene.idx,],
                            plot = FALSE, span=0.3,
                       x_lab = "Cell cycle position \u03b8", y_lab = "log2(keygene)")
    fit <- fit.l$pred.df
    options(repr.plot.width=8, repr.plot.height=6)
    print(ggplot(fit, aes(x=x,y=y)) + geom_path(size = 1, alpha = 1)+ labs(x = "θ", y = "smoothed gene expression", title = gene) + 
                scale_x_continuous(limits = c(0, 2 * pi), breaks = c(0, 
                    pi/2, pi, 3 * pi/2, 2 * pi), labels = paste0(c(0, 
                    0.5, 1, 1.5, 2), "π"), name = "θ")+
      #              scale_x_continuous(limits = c(0, 200), breaks = c(0, 
      #              50, 100, 150, 200), labels = paste0(c(0, 
      #              0.5, 1, 1.5, 2), "π"), name = "θ")+
      theme_bw(base_size = 14)+
      theme(legend.title=element_blank()))
}


## DWF4
pltsg('AT3G50660')


## OPS
pltsg('AT3G09070')


## CPD
pltsg('AT5G05690')


sce <- readRDS("./scRNA-seq/Published_Data/phloem_atlas_sce.rds")


meta <- read.csv('./scRNA-seq/Published_Data/phloem_atlas_cell_metadata.csv')


head(meta)


rownames(meta) <- meta$cell_id


nrow(meta)

sce

class: SingleCellExperiment 
dim: 17396 10204 
metadata(4): nucgenes genevar hvgs hvgs_vst
assays(4): counts logcounts vst logvst
rownames(17396): AT1G01010 AT1G01020 ... AT5G67630 AT5G67640
rowData names(10): ID chrom ... go_mitosis go_phloem_dev
colnames(10204): APL_AAACCCAAGTAACGTA-1 APL_AAACCCACAACCAACT-1 ...
  sAPL_TTCGGTCAGCGACAGT-1 sAPL_TTTACCATCCTCTGCA-1
colData names(38): Sample Barcode ... slingPseudotime_4
  slingPseudotime_5
reducedDimNames(36): MNN_logcounts MNN_logvst ... DIFFMAP_MNN_logcounts
  DIFFMAP_MNN_logvst
mainExpName: NULL
altExpNames(0):


seu <- as.Seurat(sce)

Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from MNN_logcounts_ to MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from MNN_logvst_ to MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from PC to PC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to PC_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from PC to PC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to PC_”
Warning message:
“Cannot add objects with duplicate keys (offending key: PC_), setting key to 'pca_logvst_'”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_logcounts_ to UMAP7logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_logvst_ to UMAP7logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_MNN_logcounts_ to UMAP7MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP7_MNN_logvst_ to UMAP7MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP7MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_logcounts_ to UMAP15logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_logvst_ to UMAP15logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_MNN_logcounts_ to UMAP15MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP15_MNN_logvst_ to UMAP15MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP15MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_logcounts_ to UMAP30logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_logvst_ to UMAP30logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_MNN_logcounts_ to UMAP30MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP30_MNN_logvst_ to UMAP30MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP30MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_logcounts_ to UMAP100logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_logvst_ to UMAP100logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_MNN_logcounts_ to UMAP100MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from UMAP100_MNN_logvst_ to UMAP100MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to UMAP100MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_logcounts_ to TSNE15logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_logvst_ to TSNE15logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_MNN_logcounts_ to TSNE15MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE15_MNN_logvst_ to TSNE15MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE15MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_logcounts_ to TSNE30logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_logvst_ to TSNE30logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_MNN_logcounts_ to TSNE30MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE30_MNN_logvst_ to TSNE30MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE30MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_logcounts_ to TSNE60logcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60logcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_logvst_ to TSNE60logvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60logvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_MNN_logcounts_ to TSNE60MNNlogcounts_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60MNNlogcounts_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from TSNE60_MNN_logvst_ to TSNE60MNNlogvst_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to TSNE60MNNlogvst_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_”
Warning message:
“Cannot add objects with duplicate keys (offending key: DC_), setting key to 'diffmap_logvst_'”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_”
Warning message:
“Cannot add objects with duplicate keys (offending key: DC_), setting key to 'diffmap_mnn_logcounts_'”
Warning message:
“Keys should be one or more alphanumeric characters followed by an underscore, setting key from DC to DC_”
Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to DC_”
Warning message:
“Cannot add objects with duplicate keys (offending key: DC_), setting key to 'diffmap_mnn_logvst_'”

seu

An object of class Seurat 
17396 features across 10204 samples within 1 assay 
Active assay: originalexp (17396 features, 0 variable features)
 36 dimensional reductions calculated: MNN_logcounts, MNN_logvst, PCA_logcounts, PCA_logvst, UMAP7_logcounts, UMAP7_logvst, UMAP7_MNN_logcounts, UMAP7_MNN_logvst, UMAP15_logcounts, UMAP15_logvst, UMAP15_MNN_logcounts, UMAP15_MNN_logvst, UMAP30_logcounts, UMAP30_logvst, UMAP30_MNN_logcounts, UMAP30_MNN_logvst, UMAP100_logcounts, UMAP100_logvst, UMAP100_MNN_logcounts, UMAP100_MNN_logvst, TSNE15_logcounts, TSNE15_logvst, TSNE15_MNN_logcounts, TSNE15_MNN_logvst, TSNE30_logcounts, TSNE30_logvst, TSNE30_MNN_logcounts, TSNE30_MNN_logvst, TSNE60_logcounts, TSNE60_logvst, TSNE60_MNN_logcounts, TSNE60_MNN_logvst, DIFFMAP_logcounts, DIFFMAP_logvst, DIFFMAP_MNN_logcounts, DIFFMAP_MNN_logvst


saveRDS(seu,"./tricycle/phloem_atlas_seu4.rds")


seu <- CreateSeuratObject(seu@assays$originalexp@counts, meta.data=meta)


seu <- NormalizeData(seu)
seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = nrow(seu))
all.genes <- rownames(seu)
seu <- ScaleData(seu, features = all.genes)
seu <- RunPCA(seu,features = VariableFeatures(object = seu))

Centering and scaling data matrix

PC_ 1 
Positive:  AT5G16130, AT1G27400, AT2G37190, AT3G45030, AT5G27850, AT4G29410, AT2G34480, AT3G16080, AT3G62870, AT5G22440 
	   AT1G22780, AT2G44120, AT3G11940, AT5G15200, AT1G09690, AT3G61110, AT5G20290, AT4G36130, AT2G39460, AT2G01250 
	   AT3G22230, AT5G59850, AT3G07110, AT5G60670, AT5G48760, AT1G69620, AT5G52650, AT5G39740, AT1G09590, AT5G45775 
Negative:  AT1G76180, AT1G62480, AT5G65020, AT1G65930, AT2G02130, AT1G49240, AT5G56540, AT2G19760, AT3G21770, AT3G13520 
	   AT1G10682, AT2G45820, AT1G20440, AT3G14310, AT1G05850, AT3G53420, AT5G17920, AT3G05890, AT1G12080, AT2G44790 
	   AT2G30870, AT5G11970, AT3G23820, AT3G11930, AT3G56240, AT5G28050, AT1G35720, AT2G47730, AT1G55330, AT4G26320 
PC_ 2 
Positive:  AT1G29470, AT3G58120, AT5G23860, AT2G45070, AT1G23490, AT5G14030, AT2G27030, AT1G70490, AT2G47170, AT5G15350 
	   AT3G12610, AT5G48580, AT1G78040, AT1G04680, AT1G04820, AT3G03160, AT1G10630, AT2G27130, AT2G21160, AT3G43810 
	   AT1G50010, AT5G17190, AT4G14960, AT4G24920, AT3G54810, AT2G19460, AT4G32460, AT5G03345, AT4G24190, AT3G25220 
Negative:  AT4G14010, AT2G37040, AT2G40890, AT4G34600, AT3G21240, AT5G57685, AT1G32100, AT5G64260, AT2G44300, AT5G54160 
	   AT4G34050, AT5G48930, AT3G27890, AT1G44800, AT2G36830, AT2G30860, AT1G27030, AT2G19110, AT2G40900, AT3G23090 
	   AT1G73660, AT1G18140, AT2G30490, AT3G26520, AT1G48320, AT3G10340, AT1G16350, AT1G06550, AT4G15340, AT1G48750 
PC_ 3 
Positive:  AT3G28860, AT5G19780, AT4G19410, AT1G15690, AT1G09640, AT5G60390, AT5G02500, AT3G19820, AT1G07940, AT2G31360 
	   AT4G13940, AT2G36530, AT4G38680, AT2G14890, AT4G34200, AT5G08690, AT5G14040, AT5G44340, AT3G20050, AT3G11830 
	   AT5G62890, AT5G08670, AT3G04120, AT2G34250, AT1G04430, AT1G20330, AT3G13460, AT1G24510, AT5G20890, AT5G46290 
Negative:  AT4G24960, AT1G62045, AT3G17730, AT5G18130, AT2G35585, AT5G14110, AT1G69970, AT3G49360, AT5G50290, AT3G54040 
	   AT5G02140, AT3G01940, AT1G33475, AT2G18380, AT1G18210, AT2G27740, AT1G06490, AT4G17510, AT1G65910, AT3G02677 
	   AT3G18670, AT1G63310, AT1G11570, AT1G75720, AT3G08690, AT4G05220, AT1G05610, AT1G31935, AT3G03170, AT5G01070 
PC_ 4 
Positive:  AT1G22840, AT4G32470, AT1G31812, AT1G72020, AT1G51650, AT4G21105, AT2G31490, AT5G64350, AT4G29480, AT1G15120 
	   AT3G46430, AT5G47030, AT5G53560, AT4G20150, AT1G27450, AT3G52730, AT4G16450, AT5G08060, AT3G52300, AT4G34700 
	   AT5G47890, AT4G37830, AT1G05205, AT5G59613, AT4G30010, AT1G76200, AT1G14450, AT3G10860, AT4G09030, AT5G05370 
Negative:  AT3G55760, AT3G24240, AT1G10760, AT5G14640, AT4G30020, AT2G37590, AT1G31880, AT1G06490, AT5G19090, AT1G65910 
	   AT1G55760, AT1G54030, AT2G17840, AT5G61960, AT3G02260, AT4G05430, AT3G54040, AT5G01890, AT3G49360, AT5G60200 
	   AT1G63300, AT2G18380, AT1G64390, AT4G39400, AT1G05200, AT1G05610, AT3G41768, AT1G72700, AT5G14110, AT2G37080 
PC_ 5 
Positive:  AT5G26260, AT3G16420, AT1G54000, AT2G41800, AT2G47140, AT3G53480, AT2G43610, AT3G13790, AT5G60950, AT5G10130 
	   AT3G15950, AT3G09260, AT5G02260, AT3G20370, AT1G53180, AT3G16460, AT2G04160, AT5G26280, AT2G17280, AT2G43535 
	   AT2G30340, AT1G66270, AT1G78340, AT5G35940, AT1G54010, AT5G01040, AT5G14730, AT3G29034, AT4G23590, AT2G45400 
Negative:  AT2G26730, AT4G20270, AT1G68430, AT1G60390, AT3G15680, AT5G60200, AT2G16850, AT4G22120, AT1G72150, AT2G45470 
	   AT3G52490, AT1G80690, AT5G25490, AT3G04670, AT5G57110, AT1G73850, AT3G63300, AT1G05470, AT2G18380, AT1G65910 
	   AT3G28455, AT2G01910, AT5G54660, AT4G36620, AT1G70670, AT1G56230, AT1G63300, AT2G44830, AT3G45610, AT3G59680


seu <- FindNeighbors(seu, dims = 1:10)
seu <- FindClusters(seu, resolution = 0.5)

Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 10204
Number of edges: 303357

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9067
Number of communities: 14
Elapsed time: 1 seconds


seu <- RunUMAP(seu, dims = 1:10)

10:59:20 UMAP embedding parameters a = 0.9922 b = 1.112

10:59:21 Read 10204 rows and found 10 numeric columns

10:59:21 Using Annoy for neighbor search, n_neighbors = 30

10:59:21 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

10:59:22 Writing NN index file to temp file /tmp/RtmpSzDebc/file16ca8b2cbb33de

10:59:22 Searching Annoy index using 1 thread, search_k = 3000

10:59:25 Annoy recall = 100%

10:59:27 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

10:59:30 Initializing from normalized Laplacian + noise (using RSpectra)

10:59:30 Commencing optimization for 200 epochs, with 388778 positive edges

10:59:37 Optimization finished


head(seu@meta.data)


DimPlot(seu, reduction = "umap", group.by = "annotation")


sce <- as.SingleCellExperiment(seu)

## Tricycle estimation
sce <- project_cycle_space(sce,ref.m = new.ref)
sce <- estimate_cycle_position(sce)

The number of projection genes found in the new data is 242.


sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage


## Tricycle position on umap
p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                       point.size = 4, point.alpha = 1) +
  theme_void(base_size = 14)
legend <- circle_scale_legend(text.size = 5, alpha = 1)
print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))


DimPlot(seu, reduction = "umap", group.by = "tricycleCCStage")


saveRDS(seu,"./tricycle/phloem_atlas_seu4_simplified.rds")


# Set the identity of interest to cell type
Idents(rc.integrated) <- "tricycleCCStage"
# Use SCT assay 
DefaultAssay(rc.integrated) <- "integrated"
de.list <- FindAllMarkers(rc.integrated,
                            max.cells.per.ident = 10000,
                            only.pos=T, 
                           test.use="wilcox")

Calculating cluster G2M


de.list <- de.list %>% filter(p_val_adj < 0.01 & avg_log2FC > 0.25)


table(de.list$cluster)

  G2M G1/G0     S 
  133     9   103


de.list %>% filter(cluster=="G1/G0")

	used	(Mb)	gc trigger	(Mb)	max used	(Mb)
Ncells	625275	33.4	1361470	72.8	1209572	64.6
Vcells	1159072	8.9	8388608	64.0	1802279	13.8

	Gene	Phase	Symbol	Araport11	Function
	<chr>	<chr>	<chr>	<chr>	<chr>
1	AT1G12845	G1/G0	NA	transmembrane protein;(source:Araport11)	NA
2	AT3G26960	G1/G0	NA	Pollen Ole e 1 allergen and extensin family protein;(source:Araport11)	NA
3	AT5G02260	G1/G0	ATEXP9	expansin A9;(source:Araport11)	member of Alpha-Expansin Gene Family. Naming convention from the Expansin Working Group (Kende et al, 2004. Plant Mol Bio)
4	AT1G33930	G1/G0	NA	P-loop containing nucleoside triphosphate hydrolases superfamily protein;(source:Araport11)	NA
5	AT5G01445	G1/G0	NA	hypothetical protein;(source:Araport11)	NA
6	AT5G67260	G1/G0	CYCD3;2	CYCLIN D3;(source:Araport11)	Encode CYCD3;2, a CYCD3 D-type cyclin. Important for determining cell number in developing lateral organs. Mediating cytokinin effects in apical growth and development.

	Gene	Phase	Symbol	Araport11	Function
	<chr>	<chr>	<chr>	<chr>	<chr>
1	AT1G16630	G2	NA	transmembrane protein;(source:Araport11)	NA
2	AT4G05190	G2	ATK5	kinesin 5;(source:Araport11)	ATK5 encodes a kinesin protein involved in microtubule spindle morphogenesis. It acts as a minus-end directed motor as well as a plus-end tracking protein (+TIP). Localizes to mitotic spindle midzones and regions rich in growing plus-ends within phragmoplasts.
3	AT1G20610	G2	CYCB2;3	Cyclin B2;(source:Araport11)	NA
4	AT5G55520	G2	NA	kinesin-like protein;(source:Araport11)	NA
5	AT4G01730	G2	NA	DHHC-type zinc finger family protein;(source:Araport11)	NA
6	AT1G44110	G2	CYCA1;1	Cyclin A1;(source:Araport11)	NA

	Gene	Phase	Symbol	Araport11	Function
	<chr>	<chr>	<chr>	<chr>	<chr>
1	AT5G10400	S	H3.1	Histone superfamily protein;(source:Araport11)	NA
2	AT5G22880	S	H2B	histone B2;(source:Araport11)	Encodes a histone 2B (H2B) protein. This protein can be ubiquitinated in planta, and this modification depends on the HUB1 and HUB2 E3 ubiquitin ligases.
3	AT5G59870	S	h2a.w.6	histone H2A 6;(source:Araport11)	Encodes HTA6, a histone H2A protein.
4	AT5G59690	S	NA	Histone superfamily protein;(source:Araport11)	NA
5	AT5G65360	S	H3.1	Histone superfamily protein;(source:Araport11)	NA
6	AT1G09200	S	H3.1	Histone superfamily protein;(source:Araport11)	NA

	G2	M	S
AAACCTGAGCGCCTCA_1	-0.028233902	0.02767077	-0.2789407
AAACCTGCACTTAACG_1	-0.085037611	-0.09186232	-0.3035347
AAACCTGCAGGACGTA_1	-0.128067505	-0.03061144	-0.2891731
AAACCTGGTATCACCA_1	-0.005055825	-0.07676504	-0.2121654
AAACCTGTCAGGTAAA_1	-0.082317857	-0.09259665	-0.1922454
AAACGGGCAATAGAGT_1	-0.154851348	-0.11182951	-0.3160343

Load reference¶

Top Genes contribute the most in each PCs¶

Phase signature¶

Zhang's list¶

Goldy's list¶

Check Seurat Cluster 0 and 4¶

Convert to single cell experiment¶

Tricycle inference¶

Permanent recipe for tricycle inference:¶

Protophloem SMART-seq data¶

tricycle reference¶

Phloem Atlas (Nature Plants 2022, 10k cells)¶

DE analysis¶

	gene	phase
	<chr>	<chr>
1	AT2G23320	Enriched_G2M
2	AT4G17340	Enriched_G2M
3	AT5G27420	Enriched_G2M
4	AT5G54290	Enriched_G2M
5	AT5G08240	Enriched_G2M
6	AT5G06980	Enriched_G2M

	Depleted_G2M	Enriched_G2M
AAACCTGAGCGCCTCA_1	-0.055274931	0.05566898
AAACCTGCACTTAACG_1	0.002564458	-0.10651000
AAACCTGCAGGACGTA_1	-0.072215184	-0.10086421
AAACCTGGTATCACCA_1	-0.161443430	0.01979224
AAACCTGTCAGGTAAA_1	0.028097562	-0.05966497
AAACGGGCAATAGAGT_1	0.092315593	-0.06260356

	SLX_12100_i701_i502	SLX_12100_i701_i503	SLX_12100_i701_i504	SLX_12100_i701_i505	SLX_12100_i701_i506	SLX_12100_i701_i507	SLX_12100_i701_i508	SLX_12100_i701_i517	SLX_12100_i702_i502	SLX_12100_i702_i503	...	SLX_17313_i728_i521	SLX_17313_i728_i522	SLX_17313_i729_i513	SLX_17313_i729_i515	SLX_17313_i729_i516	SLX_17313_i729_i517	SLX_17313_i729_i518	SLX_17313_i729_i520	SLX_17313_i729_i521	SLX_17313_i729_i522
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	...	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
AT1G01010	138.832	0.000	0	1.37359	0.665616	23.122100	0.0000	0	0.00000	153.1940	...	45.11770	0.00000	0	13.7543	58.2903	0.00	0.00000	0	0.0000	0
AT1G01020	277.664	146.593	0	33.42390	0.000000	0.000000	39.6054	0	0.00000	105.3210	...	8.01155	0.00000	0	0.0000	0.0000	0.00	0.00000	0	0.0000	0
AT1G01030	0.000	0.000	0	0.00000	7.321770	0.000000	0.0000	0	0.00000	0.0000	...	0.00000	0.00000	0	0.0000	0.0000	0.00	0.00000	0	0.0000	0
AT1G01040	0.000	0.000	0	0.00000	0.332834	0.379051	22.9967	0	0.00000	4.7873	...	0.00000	2.43814	0	0.0000	26.1001	0.00	0.00000	0	0.0000	0
AT1G01046	0.000	0.000	0	0.00000	0.000000	0.000000	0.0000	0	0.00000	0.0000	...	0.00000	0.00000	0	0.0000	0.0000	0.00	0.00000	0	0.0000	0
AT1G01050	138.832	146.593	0	72.80000	105.167000	137.974000	136.7020	0	1.94657	4.7873	...	6.32491	104.80900	0	103.1570	29.5801	329.67	6.51429	0	10.5481	0

	orig.ident	nCount_RNA	nFeature_RNA	ID	sample	reads_total	reads_unmapped	reads_unique	reads_multi	percent_unique	percent_multi	percent_mapped	lib_size_avr	Cluster_all_cells	Cluster_PSE	Pseudotime_PSE
	<fct>	<dbl>	<int>	<chr>	<chr>	<int>	<int>	<int>	<int>	<dbl>	<dbl>	<dbl>	<int>	<int>	<int>	<dbl>
SLX_12100_i701_i502	SLX	734281.9	2860	SLX_12100_i701_i502	P1	847809	803895	38147	5767	4.499480	0.6802240	5.179704	678	NA	NA	NA
SLX_12100_i701_i503	SLX	913422.6	3518	SLX_12100_i701_i503	P1	1431659	1399319	27698	4642	1.934679	0.3242392	2.258918	678	NA	NA	NA
SLX_12100_i701_i504	SLX	827476.9	2904	SLX_12100_i701_i504	P1	313281	300049	9440	3792	3.013269	1.2104149	4.223684	678	NA	NA	NA
SLX_12100_i701_i505	SLX	600775.5	10725	SLX_12100_i701_i505	P1	3328241	992332	1515355	820554	45.530206	24.6542844	70.184491	678	10	4	1.3333992
SLX_12100_i701_i506	SLX	922948.5	11363	SLX_12100_i701_i506	P1	7051643	1901022	3201903	1948718	45.406482	27.6349498	73.041432	678	14	4	0.8876946
SLX_12100_i701_i507	SLX	894200.0	10634	SLX_12100_i701_i507	P1	6906310	1881499	2873579	2151232	41.608022	31.1487900	72.756812	678	7	NA	NA

A data.frame: 6 x 9
	cell_id	sample	barcode	total_counts	detected_genes	cluster	annotation	UMAP1	UMAP2
	<chr>	<chr>	<chr>	<int>	<int>	<int>	<chr>	<dbl>	<dbl>
1	APL_AAACCCAAGTAACGTA-1	APL	AAACCCAAGTAACGTA-1	35600	5210	9	early/undetermined	2.2210797	-3.4197492
2	APL_AAACCCACAACCAACT-1	APL	AAACCCACAACCAACT-1	6616	3058	2	PSE	6.6928352	5.7139691
3	APL_AAACCCACAGTTGAAA-1	APL	AAACCCACAGTTGAAA-1	12781	3669	1	MSE	0.9346761	0.6447128
4	APL_AAACCCAGTGCCTTTC-1	APL	AAACCCAGTGCCTTTC-1	7138	3100	4	PPP	-5.9457399	1.5383996
5	APL_AAACCCAGTTTCGACA-1	APL	AAACCCAGTTTCGACA-1	7364	2529	6	PSE	5.8819741	8.5710071
6	APL_AAACCCATCTGCGAGC-1	APL	AAACCCATCTGCGAGC-1	7114	2869	3	CC	-0.7623315	-0.4224323

A data.frame: 6 x 14
	orig.ident	nCount_RNA	nFeature_RNA	cell_id	sample	barcode	total_counts	detected_genes	cluster	annotation	UMAP1	UMAP2	RNA_snn_res.0.5	seurat_clusters
	<fct>	<dbl>	<int>	<chr>	<chr>	<chr>	<int>	<int>	<int>	<chr>	<dbl>	<dbl>	<fct>	<fct>
APL_AAACCCAAGTAACGTA-1	APL	35383	5176	APL_AAACCCAAGTAACGTA-1	APL	AAACCCAAGTAACGTA-1	35600	5210	9	early/undetermined	2.2210797	-3.4197492	8	8
APL_AAACCCACAACCAACT-1	APL	6456	3020	APL_AAACCCACAACCAACT-1	APL	AAACCCACAACCAACT-1	6616	3058	2	PSE	6.6928352	5.7139691	9	9
APL_AAACCCACAGTTGAAA-1	APL	12731	3650	APL_AAACCCACAGTTGAAA-1	APL	AAACCCACAGTTGAAA-1	12781	3669	1	MSE	0.9346761	0.6447128	8	8
APL_AAACCCAGTGCCTTTC-1	APL	6954	3072	APL_AAACCCAGTGCCTTTC-1	APL	AAACCCAGTGCCTTTC-1	7138	3100	4	PPP	-5.9457399	1.5383996	4	4
APL_AAACCCAGTTTCGACA-1	APL	7276	2483	APL_AAACCCAGTTTCGACA-1	APL	AAACCCAGTTTCGACA-1	7364	2529	6	PSE	5.8819741	8.5710071	11	11
APL_AAACCCATCTGCGAGC-1	APL	7045	2849	APL_AAACCCATCTGCGAGC-1	APL	AAACCCATCTGCGAGC-1	7114	2869	3	CC	-0.7623315	-0.4224323	1	1

A data.frame: 9 x 7
	p_val	avg_log2FC	pct.1	pct.2	p_val_adj	cluster	gene
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<fct>	<chr>
AT5G02820	1.728994e-120	2.7087332	0.138	0.378	4.236034e-118	G1/G0	AT5G02820
AT1G04710	2.012929e-90	2.8375044	0.120	0.281	4.931675e-88	G1/G0	AT1G04710
AT2G03670	4.696481e-50	4.8883449	0.138	0.325	1.150638e-47	G1/G0	AT2G03670
AT3G13060	2.300513e-25	1.9837473	0.133	0.147	5.636256e-23	G1/G0	AT3G13060
AT2G42190	1.941702e-21	0.8564682	0.126	0.202	4.757170e-19	G1/G0	AT2G42190
AT4G14980	3.124803e-08	2.9606000	0.111	0.129	7.655767e-06	G1/G0	AT4G14980
AT1G05520	1.283976e-07	1.8709291	0.183	0.239	3.145741e-05	G1/G0	AT1G05520
AT1G128001	1.191791e-06	0.7292746	0.141	0.264	2.919888e-04	G1/G0	AT1G12800
AT5G63120	2.330971e-06	4.8604293	0.110	0.179	5.710880e-04	G1/G0	AT5G63120